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1. Introduction 

Semiclassical laser theory, which neglects the quantum fluctuations of the electromagnetic field, 
is widely used to describe and simulate lasers dO. In principle, it correctly describes the laser 
thresholds and frequencies, the spatial pattern of the lasing modes, and the laser output power, 
including all classical non-linear effects, such as spatial hole-burning, gain saturation, and mode 
and phase locking. Essentially the theory describes Maxwell's equations in an open cavity, cou- 
pled to the non-linear polarization of the gain medium. The gain polarization can be described 
using either a classical non-linear oscillator model O , or a quantum-mechanical model of N 
atomic levels in which the polarization and level populations obey the equations of motion of 
the quantum density matrix. The simplest version of the theory, used widely in textbooks, is the 
two-level Maxwell-Bloch (MB) model [1|; however, most design and characterization simula- 
tions of lasers use models with N = 3 or more levels. In addition, most theoretical solutions for 
the semiclassical laser equations employ a large number of simplifying assumptions in order 
to make them analytically tractable, most notably neglecting the openness of the cavity and/or 
treating only simple one-dimensional (ID) or ring cavities, as well as approximating the non- 
linear interactions to cubic order. The results are typically not useful for quantitative modeling. 
Until recently, the only useful way to obtain quantitative results for non-trivial laser structures 
was to integrate the semiclassical laser equations in space and time. For novel and interest- 
ing new laser structures with non-trivial 2D and 3D cavity geometries, such simulations are at 
the limits of computational feasibility, making it difficult to study a large parameter space or 
ensemble of designs. 

In the past five years a new approach to finding the stationary solutions of the semiclassical 
laser equations has been developed, known as Steady state Ab initio Laser Theory (SALT) (61 [71 
migl. SALT treats the openness of the cavity exactly, and the multimode non-linear interactions 
to infinite order (with two approximations, to be discussed below). It is applicable to cavities 
of arbitrary complexity in 2D and 3D, although we will discuss here only the scalar wave 
equation coupled to a gain medium. Importantly, this approach eliminates the need to perform 
a time integration to steady-state, dramatically reducing the computational effort and allowing 
one to study cavities with high spatial complexity, such as 2D random lasers |[8l[T0l|, photonic 
crystal lasers (HI [111 and chaotic disk microlasers (T3l[T4ll . 

SALT was originally formulated to solve the steady- state two-level Maxwell-Bloch equa- 
tions (51 [71 in the standard slowly- varying envelope approximation (SVEA), and an iterative 
algorithm was developed |7 1 to solve the resulting SALT equations. Subsequently it was real- 
ized that the SVEA afforded no advantage in numerical solutions of the SALT equations, so 
this approximation was dropped, leading to slightly different SALT equations (l4l[T5l, which 
were used in all subsequent work [Si Qi il2j| . Recently an important generalization of the SALT 
solution algorithm was developed, improving its performance for laser cavities with complex 
spatial index variation and inhomogeneous pumping f9\ . In its present form, SALT only em- 
ploys two approximations, the stationary inversion approximation (SIA) and the rotating wave 



approximation (RWA), which are both well satisfied for most lasers of interest. In Ref. lITSl , 
the results of SALT were compared to the full time-dependent solution of the MB equations for 
a simple ID cavity in the multimode regime, well above threshold. Excellent agreement was 
found in the parameter regime for which the SIA holds. This was, to our knowledge, the first 
demonstration of a frequency-domain method which agrees with exact time-domain methods 
for above-threshold multimode lasing. 

Previous applications of SALT have focused on two-level gain media. In order for SALT to 
be a useful modeling tool, it is necessary to demonstrate that it can be applied to A/^-level lasers. 
In the current work, we show analytically that the steady- state equations for an A/^-level laser 
can be reduced to those for an effective two-level system, and hence solved using the efficient 
SALT algorithm with essentially the same degree of computational effort. We also explore 
how this effective two-level system differs from the ordinary two-level laser. Next, we present 
a numerical comparison between the results of SALT calculations and exact A/^-level finite- 
difference time-domain (FDTD) calculations, for the same simple ID laser studied in ifTSl . as 
well as for a ID random laser. We note that a similar comparison between SALT and FDTD 
has been performed for a four-level high-Q single-mode photonic crystal laser in Ref. [12J, with 
good agreement found. Here, we test SALT'S accuracy in treating the more challenging case of 
multimode lasing in low-Q and random cavities. 

2. Effective two-level systems 

2.1. Four level system analysis 

We illustrate the approach using the semi-classical laser equations fTl for the four-level atomic 



gain medium shown in Fig.[TJ 

4;rP+ = c^V^E+-ec(r)E+ (1) 

P+ = -(/co, + n)P+ + ^E+(p22-pii) (2) 

in 

P33 = ^(P00-P33)-r23P33 (3) 

P22 = r23P33-ri2P22-^E+ •((?+)*-?+) (4) 

Pll = ri2p22-r01Pll + ^E+ •((?+)*-?+) (5) 

Poo = r0lPll-^(P00-p33). (6) 



The RWA has already been made, and we have assumed a lasing structure with one or two direc- 
tions of translational symmetry, so that TM and TE polarizations are conserved and Maxwell's 
equations reduce to a scalar Helmholtz equation |[T6t . E+ and P+ are the positive frequency 
components of the scalar electric and polarization fields respectively, pu is the population den- 
sity of level I /) , (Oa is the frequency of the gain center, yx_ is the gain width (polarization dephas- 
ing rate), g is the dipole matrix element, 3^ is the pump rate, and jij is the decay rate from level 
I j) to level |/). The four levels are labelled from — 3 in order of increasing energy (Fig. [T]). 
The polarization equation ^ is obtained from the four-level density matrix equation of motion, 
assuming that only the level 2^1 transition will be inverted and lase. Often in FDTD calcu- 
lations a real classical oscillating dipole equation is used to describe the polarization I12iil ; in 
Appendix A, we show that this yields essentially the same results, with an appropriate identifi- 
cation of parameters. In the rate equations ([S])-®, the pump is coherently acting between levels 
|0) and |3). 

The polarization equation Q incorporates the inversion D{x^t)= p22 (•^, ~ Pi i (-^^ ' similar 
to the polarization equation for a two-level laser. By assuming that the non-lasing populations 
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Fig. 1. Schematic of a four-level gain medium. 



are stationary, poo = P33 = 0, we can show that D{x,t) obeys 

Z) = -T^(D-D;,)-|e+ •((?+)*-?+), (7) 

which is precisely the form of the inversion equation for the two-level medium JT] [181 • The 
parameters Dq and serve as an effective equilibrium inversion and inversion relaxation rate 
respectively, and are given by llT4l : 

where S = {joi — 7i2)/7i2 and n = Y^iPu is the total density of gain atoms. Eq. © for the 
inversion in the absence of laser emission (which here acts as the effective pump parameter) has 
been discussed by Siegman ||3J, while Eq. ([5]) for the effective relaxation rate has been derived 
for a special case by Khanin |[T9l . These expressions have not been used previously to solve 
the four-level lasing equations in terms of the two-level solutions, as we do here. If we use an 
incoherent pump, the ^(poo — P33) terms in ^ and © would be replaced by ^poo', the effect 
of this would be the removal of the factor of 2 preceding the ratio 701 / 723 in the denominators 
of ([8]) and (|9]). For typical laser systems, each denominator is dominated by the 701 / ^ term, so 
the difference between coherently and incoherently pumping the system is negligible 1 17 |. 

Given the parameters {701 , 712, 723, ^} describing the four-level medium of Fig. [2l we can 
calculate the effective pump and relaxation rates, and feed those (together with the cavity dielec- 
tric function, lasing transition frequency cOa, and gain linewidth 7^) into the SALT algorithm. 
That will yield the steady-state lasing properties of this four-level laser. 

2.2. Arbitrary number of levels 

A general A/^-level laser can be treated via an analogous procedure. Suppose we have a gain 
medium with an arbitrary number of levels, N. We assume that there is only a single lasing 
transition, between two levels which we denote by \u) and |/). The population of each of the 
N — 2 non-lasing levels obeys 



Pa = Y.yiJpJj -HyjiPii + yiuPuu + YuPu - {Yui + rii)pii. 

j j 



(10) 



where the sums are taken over all of the non-lasing states and Jij is the rate at which atoms 
transition from state |j) to |/) which is either a decay or pump rate depending on the relative 
energies of the states. In this way, we can incorporate decay and pump processes between any 
levels. We again assume that pa = for all non-lasing levels. Then 

52 [(^i + + Jld^ij - Jij] Pjj = YiuPuu + YilPlh (11) 

j 

where st = ^jYji and 5ij is the Kronecker delta. The term in brackets on the left hand side 
corresponds to an {N — 2)x {N—2) matrix, which we denote as R. Upon inverting Eq. (fTTl) and 
substituting it into the equations of motion for the lasing levels, we obtain 

, _BiTu-BuTi , Bu^Bi n 

^1- 7;+7] ' ^^-j-pFij^ ^^2^ 



where 

Tu/i = ^+"L[R~%jyj,uii (13) 

ij 

Bu = -su^Y.^yui-yii)[R~%jyjm (i4) 

ij 

Bi = si^J^{rui-rii)[R-']ijrji- (i5) 

ij 

The details of this calculation are given in Appendix B. 
3. Physical Limits of Interest 

Returning to the typical four-level case, we take note of two important physical regimes. The 
first, the linear regime, is 723 ~ 7oi ^ 712 ^ for which one recovers the expected behavior 
that the equilibrium inversion increases linearly with the pump and that is a constant: 

7y - ^Yn. (16) 
D'q « —n. ill) 

7l2 

In this case, varying the equilibrium inversion and the pump strength are essentially equivalent. 

The second regime of interest, the non-linear regime, is 7^3 ~ 701 ^ 712 ~ i-e. when the 
slow decay rate between the lasing levels is on the same order as the pump rate. In this regime, 
7| increases with increasing pump and D'^ saturates with increasing pump: 

0^ « 2(712 + ^), (18) 
1 / ^\ 

D'o « : ^ — «■ (19) 



" ^ ^V7.. 

This regime is also interesting from the viewpoint of SALT. As 7y increases with a laser 
could satisfy the inequality 7y <C 7_l near threshold, leading to stationary inversion and an accu- 
rate solution via SALT, but fail to satisfy the inequality as the pump becomes stronger, leading 
to a decrease in the accuracy of SALT. 

For a system with an arbitrary number of levels, the first regime always occurs at sufficiently 
small pump values; the second regime is obtainable if electrons in the upper lasing level are 
relatively long-lived compared to electrons in other levels. 



4. Brief Summary of SALT 

For completeness we briefly outline SALT. The and fields are assumed to obey a multi- 
mode ansatz 

M 

£+(r,0= £l'^(?)e-''»'", 

(20) 

M ^ ^ 

11=1 

where the indices /i = 1 , 2, • • • , M label the different lasing modes, and the field and polarization 
are now explicitly scalar quantities. The total number of modes, M, is not given, but increases 
in unit steps from zero as we increase the pump strength Dq. The values of Dq at which each 
step occurs are the (interacting) modal thresholds, to be determined self-consistently from the 
theory. The real numbers (0^ are the lasing frequencies of the modes (henceforth c = 1), which 
will also be determined self-consistently. 

We insert the ansatz (l20l) into the two-level laser equations, and employ the stationary inver- 
sion approximation (SIA) D = 0. The result is a set of coupled nonlinear differential equations, 
which are the fundamental equations of SALT ||9l : 



Dir)=Do{7) 



M 1 

i + £rv|>Pv(r)p 

v=l 



^m(^)=0. (21) 

(22) 



^ and D are now dimensionless, measured in their natural units Ec = fi-sJy\{Yl/ (2g) and 
J^c = ^7±/(47rg^), and Fy = y\/{y\ + (cOy — (OaY) is the Lorentzian gain curve evaluated 
at frequency (Oy. Note that these equations are time-independent; Eq. ([2T]) is a stationary wave 
equation for the electric field mode , with an effective dielectric function consisting of both 
the "passive" contribution Ecij) and an "active" contribution from the gain medium. The latter 
is frequency-dependent, and has both a real part and a negative (amplifying) imaginary part. 
It also includes infinite-order nonlinear "hole-burning" modal interactions, seen in the |^yp 
dependence of (1221) . In addition, we make the key requirement that must be purely out- 
going outside the cavity; it is this condition that makes the problem non-Hermitian. It is worth 
noting that the SIA is not needed until at least two modes are above threshold, so (l22l) is exact 
for single-mode lasing up to and including the second threshold (aside from the well-obeyed 
RWA). 

These equations are solved efficiently by projecting them onto a complete biorthogonal set of 
purely outgoing states with external wavevectors, Z:^, equal to the lasing frequencies. We refer 
to these states as the threshold constant flux (TCF) states, because one member of the basis set is 
always equal to the (non-interacting) threshold lasing mode, leading to very rapid convergence 
of the basis expansion above threshold (91. The major computational effort in solving the SALT 
equations by this approach is in calculating the (linear) TCF states. The SALT solutions are 
obtained for successive values of the pump increasing from the first threshold; at each step, 
the coefficients of the lasing modes in the TCF basis are obtained using a standard non-linear 
solver, with the coefficients from the previous step as an initial guess (which is never far from 
the correct solution). Unlike FDTD, in the current version of SALT one cannot simply directly 
solve at a fixed pump value, well above threshold. However, even with this limitation, SALT is 
much more efficient, and provides substantial physical insight, as we will discuss below. 



5. Numerical comparison 

To perform a well-controlled comparison between SALT and the four-level laser equations ©- 
([6]), as well as A/^-level generalizations, we studied ID microcavity lasers for which the FDTD 
calculations are tractable and fast enough to generate extensive steady- state data. We first con- 
sider the same simple edge emitting uniform-index laser treated in Refs. (61 Ul El, with a 
perfect mirror at the origin, active region of length L terminating abruptly in air (see schematic. 
Fig. O. The simulations were carried out using standard FDTD for the electromagnetic field, 
and Crank-Nicholson discretization for the polarization and rate equations based on the method 
of Bidegaray ll20t (in which the polarization and inversion are spatially aligned with the electric 
field but updated at the same time steps as the magnetic field). The reported modal intensities 
are calculated by Fourier transforming the electric field at the cavity boundary after the simu- 
lation has reached steady state (see ^Hfor a discussion of the steady-state criterion). The lasing 
transition frequency cOa is chosen so that nokaL = 60, corresponding to roughly ten wavelengths 
of radiation within the cavity. Physical quantities are reported in terms of their natural scales, 
Dc = hy±/ (47rg^) and Ec = (^/2g) y^TLT^- addition, we take c = h = l and measure rates in 

dimensionless units, i.e. /meas = 7rQa\L/c. We note that the parameters chosen accurately reflect 
those of real microcavities at optical frequencies [12 1 ; the complete set of simulation parameters 
is given in Appendix D. 
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Equilibrium Inversion Dq/Dc 

Fig. 2. Modal intensities as functions of the normalized equilibrium inversion Dq/Dc (ef- 
fective pump) in a ID microcavity edge emitting laser (schematic inset). The cavity is 
bounded on one side by a perfect mirror and on the other side by air, and has uniform 
refractive index n = 1.5. Solid lines are results obtained by the time-independent SALT 
method; open circles are results of FDTD simulations with a coherently pumped four-level 
medium (Fig.[T]); solid triangles are results of FDTD simulations with a coherently pumped 
six-level medium with a lasing transition between |3) and |1). Simulation parameters are 
given in Appendix D. Both the four-level and six-level media are chosen to satisfy SIA. 
The dephasing rate is yj_ = 4.0. The four-level system is in the linear regime described by 
Eq. ([T6l)-([T7l). The six-level system is calculated using the formula in the appendix B, but 
is in the non-linear regime described by Eq. ([T8]|-([T9]|. The spectra at Dq/Dc = 0.488, and 
the gain curve, are shown in the upper left inset. 

As shown in Fig. [21 we find close agreement between SALT calculations and FDTD simula- 
tions. At a representative pump strength Dq = 0.488/)^, the mode intensities produced by SALT 
differ from those of the four- and six-level FDTD simulations by ^ 1%, while the frequencies 



differ by < 0.1%. The difference in mode frequencies between SALT and FDTD also exists at 
the first lasing threshold, for which an analytical value can be calculated. There, we find that 
the FDTD simulation has a 0.2% error in the first mode frequency, while SALT has a 0.08% 
error; this error arises from the spatial discretization of the cavity employed in both approaches 
II2TII . It is worth emphasizing that SALT treats the non-linearity to infinite order; in the earlier 
work on the Max well-B loch model ifTSll it was shown that for this same cavity the common 
cubic approximation for the non-linearity fails both quantitatively and qualitatively. 

These results demonstrate that so long as the system satisfies SIA, the mapping between 
systems with an arbitrary number of levels to an effective two-level system is nearly exact, 
and SALT is able to very accurately determine the steady state properties of the cavity. If two 
cavities, each with an arbitrary number of levels, have the same effective parameters D'^ and , 
and otherwise have the same polarization relaxation rate and atomic transition frequency, the 
cavities are equivalent from the electromagnetic point of view, and will have identical lasing 
properties. 

The six-level simulations shown in Fig. [21 occupy the non-linear parameter regime of 
Eq. ([T8l)-([T9b, i.e. 7y is a linear function and Z)q a non-linear function of 3^. However, the 
unsealed modal intensity leaving the cavity is still, to leading order, linear in 3^. This can be 
seen by rearranging (1221) . inserting the expressions for 'j^ and Z)q, and noting that at the end of 
the cavity the inversion is roughly independent of the pump strength. This result is discussed 
further in Appendix C. 



Fig. 3. Breakdown of the equivalence between SALT and FDTD when SIA is not valid 
is shown here in two different ways. Here, modal intensities as a function of the normal- 
ized equilibrium inversion D'^/Dc (effective pump) are shown for a ID microcavity edge 
emitting laser with 7^ = 4.0 and n= 1.5. Solid lines again represent results obtained from 
SALT, while open circles represent FDTD simulations of a simple four-level system with 
Tjl =0.1. Triangles represent FDTD simulations of a six-level system in the non-linear pa- 
rameter regime in which 7^ ~ 0.001 for Z)q < 0.1, and thus satisfying SIA, but 7|| ~ 0.01 
for D'q > 0.45, and consequently no longer satisfying SIA. 

The mapping between the A/^-level laser and two-level SALT breaks down at large pump 
strengths, when the condition 7y <C 7_l is violated due to the increase of 'j^ with In Ref. ifTSll . 
following an argument by Haken (HI, it was demonstrated that violating this condition causes 
the SIA for the two-level model to break down. This effect can be seen in the four- level laser 




data in Fig. [51 where = 0.1, = 4.0 and accuracy is already lost for the third lasing mode. 
For the six-level data of Fig. [51 which is in the non-linear parameter regime, SIA is satisfied 
and SALT agrees with the FDTD simulations for small values of the normalized equilibrium 
inversion; for larger values of Z)q, the SALT and FDTD results begin to diverge. 

Finally, to demonstrate that the mapping to an effective two-level model works equally well 
for a complex laser cavity. Fig. [H shows a comparison between SALT and FDTD simulations 
for a four-level gain medium in a ID random dielectric structure. A number of studies have 
been published on random lasers using such simulations |[22l [251 ; SALT provides a much more 
efficient method for such studies, which often require generating a statistical ensemble of lasers. 
Here, the passive cavity dielectric function contains ~ 31 layers, alternating randomly between 
regions with refractive indices = 1.25 and n2 = l. Each random layer was generated accord- 
ing to the formula Ji, 2 = (Ji,2)(l + ^C) where {di) = (l/3)(L/30) and (J2) = (2/3)(L/30) are 
the average thicknesses of the layers, t] = 0.9 represents the degree of randomness of the cavity, 
and G [ — 1 , 1] is a randomly generated number. The gain medium was added uniformly to the 
entire cavity, and the coherent pump was likewise uniform. The transition frequency was cho- 
sen such that nikaL = 120, corresponding to roughly 20 wavelengths inside of the cavity. We 
find only small discrepancies between the SALT and FDTD results, with ^1.1% difference in 
the modal intensities. These differences did not vary significantly between different realizations 
of the random laser. 
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Fig. 4. SALT and FDTD results for a ID random laser. Modal intensities are plotted against 
the normalized equilibrium inversion Dq/Dc (effective pump). Solid lines represent SALT 
results, and circles represent FDTD simulations for a four-level system with = 0.001. 
The refractive index distribution of the edge emitting random laser is described in the text. 
The gain medium has Yj_ = 4.0 and is in the regime described by Eq. ([T6l-([T7l. Left inset: 
log-log plot of the indicated region where three modes turn on in close proximity. Right 
inset: schematic of the cavity structure. 




6. Computational efficiency of SALT for N-level systems 

In this section we present a set of benchmarks comparing the computational efficiency of SALT 
to FDTD. SALT calculations enjoy three main advantages over FDTD simulations of the semi- 
classical laser equations. First and foremost, SALT directly finds the steady-state solutions, so 
no time integration is involved, which substantially decreases computational effort. Second, 
SALT unambiguously determines how many modes are lasing at a given pump, whereas it 



can be difficult to determine, especially for multimode lasing, when an FDTD simulation has 
reached the steady-state with all modes that will lase "on". Third, within SALT, with minimal 
additional computational effort, it is possible to monitor modes which are below threshold via 
a modified threshold matrix O, and hence to ascertain if more modes are likely to turn on in 
some interval of pump. It is important to note that the implementation of SALT used in this 
study has also not yet been fully optimized. For instance, the implementation of SALT used in 
this study requires calculating the entire lasing intensity and frequency spectrum starting from 
the first lasing threshold. However, this is not necessary and it is possible to implement SALT 
to take an initial guess of the number of modes and their relative intensities and then allow the 
algorithm to flow to the correct solution as the SALT algorithm has been demonstrated to be 
rather robust |[T4ll . Thus, while one might assume from Fig. [5] that there would be a crossover 
pump value, high above threshold, where it would be more efficient to calculate the steady-state 
solutions using FDTD, this is not necessarily the case. 

SALT does have one disadvantage that FDTD does not, assuming the cavity is not in the 
chaotic regime. The convergence time in FDTD is determined by the longest time scale in the 
problem, which is the greater of beating period between consecutive modes and the relaxation 
oscillation time. These time scales are relatively independent of the number of modes lasing in 
the cavity, so the efficiency is largely independent of the pump in the multimode regime. This 
is not the case for SALT, as the computational time increases as where N is the number of 
lasing modes. However, as we see in Fig. [5j even a non-optimal implementation of SALT is 
substantially more efficient than FDTD even when calculating the steady state of a single pump 
value. 
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Fig. 5. Comparison of SALT and FDTD run-times. Modal intensities are shown as a func- 
tion of the run-time for SALT (squares) and four-level FDTD simulations (circles), using 
the parameters of Fig. [51 FDTD simulations that have not begun to lase are marked as 
crosses. Plot (a) shows data for D^jDc = 0.071, just above the first lasing threshold. SALT 
determined the steady-state single modal intensity in under three minutes, while the FDTD 
required ^ 5000 minutes to reach steady state. Plot (b) shows data for Dq/Dc = 0.486, 
well above the third lasing threshold. SALT calculated all data up to and including this 
pump value in under 90 minutes, whereas FDTD required > 500 minutes for the first two 
modes to reach steady- state, with the third mode intensity (green circles) still fluctuating 
after 5000 minutes (not shown). 



Calculating full modal intensity/frequency curves as a function of the pump strength, such 
as in Fig. [21 is generally much more efficient using SALT. For example, in order to generate the 
curves seen in Fig. [21 SALT ran for a little under 2 processor hours. To generate all of the FDTD 
data for the four level simulations took 267 processor days. If one is attempting to explore a 



large parameter space of designs or system parameters, SALT may make studies feasible which 
are simply impractical using FDTD, particularly in more realistic 2D and 3D structures. 

As mentioned before, the bulk of the computational effort required for the SALT algorithm, 
especially in higher dimensions, is in solving for the TCF states. While the difficulty of solving 
for the TCF states does scale with the dimensionality of the system, one only need solve the 
associated generalized eigenvalue problem ~ 100 (even in higher dimensions) to have a suf- 
ficiently complete basis for all pump values, whereas in FDTD one needs to solve an 0{n^) 
problem at each of many thousands, if not millions, of time steps. Furthermore, it is likely 
that switching to a finite element method in higher dimensions is not only possible, but likely 
to result in increased computational efficiency. Once one has a TCF basis library, using the 
SALT algorithm to iterate above threshold does not directly scale with the dimensionality of 
the system. Finally, while current implementations of SALT assume the electric field is per- 
pendicular to the direction of wave propagation, the vectorial generalization of SALT is under 
investigation. 

7. Conclusion 

We have found that using stationarity conditions on the non-lasing level populations, the rate 
equations for an A/^-level laser can be mapped to an effective two-level model, for which the 
steady- state multimode lasing properties are efficiently solvable using Steady- state Ab initio 
Laser Theory (SALT). Using this mapping, we found excellent agreement between SALT and 
FDTD simulations for the modal frequencies, thresholds and above-threshold intensities, in 
the expected domain of validity. SALT is typically several orders of magnitude more efficient 
computationally than time domain solution of the laser-rate equations, assuming only steady- 
state properties are needed. 
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A. Classical polarization in the laser-rate equations 

Throughout this paper, we have used the density matrix equations of motion for the polarization 
field. However, much of the literature uses the classical oscillating dipole equation, in which the 
gain atoms are assumed to be dipoles undergoing harmonic oscillations and the quantum density 
matrix is neglected. This appendix briefly demonstrates that the two models are equivalent, and 
derives the relevant parameter redefinitions. A more thorough discussion has been given by 
Boyd El. 

The polarization equation used here. 




'uu 



pii) 



(23) 



is derived from the density matrix equation 



Pul = {i(Oa + r±) Pul - jgEip, 



'uu 



Pll) 



(24) 



where \u) and |/) are the upper and lower lasing states, pij is the density matrix element ij and 
gul = 8lu = 8 is the coupling constant of the lasing states to the electric field, which allows for 
the definition = gpi^i. Alternatively, following the notation of Boyd [.24] . we could consider 



the equation of motion for M = gp^/ + c.c. = P+ + P , which is the expectation value of the 
dipole moment induced by the appHed field, i.e. the classical oscillating dipole field. Thus 



using (l24l) . and 



M = gpui + c.c. = {icoa + r±) Pul - jgED + c.c, (25) 

n 



M= (-0)^2 + 2/0)^7^ + ri)gp,/--^g2^Z) + c.c. (26) 



This can be rewritten as 

M + 27_lM + co^M = - y\M - ^^^£Z), (27) 

n 

where D = puu — Pii is the inversion density of the lasing states. For ((i^^ we can discard 
the term 7^M, resulting in the traditional form of the classical oscillating dipole polarization 
field equation. This gives the classical coupling constant a = 2cOag^/h |3|. 

A similar analysis can used to show that the inversion equation can be rewritten as 

2E 

I) = -]^|(Z)-Z)o) + — M, (28) 

nCOa 

which demonstrates how the change in the inversion is dependent upon the classical polarization 
field. Therefore, so long as (0^ 7^, a condition that is usually satisfied, these two formulations 
of the polarization dynamics are equivalent. 

B. Effective two-level parameters from N-level rate equations 

This appendix derives the effective equilibrium inversion and relaxation rate for a gain medium 
with an arbitrary number of levels. We allow decays between any two levels, even if they are 
not adjacent. We assume there is a single lasing transition, between levels \u) and |/), which 
need not be adjacent. The rate equation for an arbitrary non-lasing level in the system is 

Pa = L yijpjj ~ L yjipii + yi^p^^ + yiipii ~ ^y^i + ^^^■)/^« ' 

j j 

where the summations are taken over all non-lasing levels. Here we do not distinguish between 
decay rates and pumping rates; jij is simply interpreted as the rate at which level I7) transitions 
into level |/), regardless of the energies of those states. If the populations of all the non-lasing 
transitions are stationary, i.e. pa = 0, then we can rewrite ([29b as 

X^^ijPjj = JiuPuu^JilPlh (30) 
j 

Rij = {si^rui^rii)Sij-rij- (3i) 

Here, Si = ^jYji and 5/^ is the Kronecker delta. Inverting (l3Qb gives 

Pa = Y^i^'^iYjuPuu + YjiPii)' (32) 

j 

Hence, we can express the total number density of gain atoms as 

^ = 1^P«+Pmm + P// (33) 

i 

= TuPuu + TiPii (34) 



where 



Tu = l^Y.[R-']ijrju. (35) 

ij 

Ti = l^Y^lR-'hjYji. (36) 

Noting that D = puu — piu we can write the populations of the lasing states as 

_ n^TiD 

J-l^ -l-u 

n-TuD 

Pll = rj. , rj. • (38) 

^/ + ^M 

From the equations of motion for the lasing levels, we have the inversion equation 

D = Puu- Pll = Y^irui - rii)pn - Supuu + sipii - • ((P+)* - P+) , (39) 
where Su = ^jjju + Jiu and si is defined similarly. Inserting ([32l) into this equation gives 

b = BuPu^u^BiPi^i - • ((P+)* - P+) , (40) 

where 

Bu = -su^Y.{rui-rii)[R~%rju, m 

ij 

Bi = si^Y^irui-Yiim-'hjYji- (42) 

ij 

Plugging in ([37l) and ([38l) now yields the inversion equation in the desired form d?]), with 

D'o = -^f^n. (44) 
C. Inversion as a function of the pump 

In this appendix we discuss the surprising result that the unsealed modal intensities of the 
six-level simulations discussed in Fig. [2] as a function of the pump are approximately linear, 
as shown in Fig. [6l even though these six-level simulations are in the non-linear parameter 
regime. In simple treatments of lasers O the inversion is assumed to be clamped after the first 
lasing threshold. However, a more detailed analysis lead to the inclusion of spatial hole-burning 
effects which forces the inversion in the cavity to change beyond the first lasing threshold in 
a non-linear manner. This result then, that the unsealed modal intensity is linear in the pump 
rate, can be understood from the observation that at the end of the cavity, r = L, all of the lasing 
modes have a maximum in their fields, and thus the inversion is effectively clamped at this point 
beyond the first lasing threshold, which can be seen in Fig.[6](b). 

To understand this formally, we begin by rewriting (l22l) and removing the scaling factors Ec 
and Dc gives 



i;r,|*,«p = ^(^j^-i)^ (45) 



Substituting in ([TSl) and ([T9l) , which are valid for this simulation, gives 



■'f;rv|^v(F)p = 



v=l 



W) 



712. 



(46) 



The inversion Z) is a function of both position and the pump. However, for r = L corresponding 
to the cavity edge, D should be mostly independent of the pump, as at this location every mode 
is at its maximum intensity and the effect of spatial hole-burning is most pronounced. The 
FDTD simulation results, shown in Fig.[6l demonstrate that D indeed varies very weakly with 
^ at the cavity edge. 
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Fig. 6. (a) Unsealed modal intensity of the six-level simulations from Fig. [2] as a function 
of the pump. A cross sectional area of Im^ is assumed to calculate the power, (b) Inversion 
as a function of the pump at the cavity boundary. Dashed lines in plots a and b correspond 
to the pump values shown in plot c. (c) Inversion as a function of position within the cavity 
for three different pump values, cyan corresponds with ^ = 3.75 x \0^s~^ , magenta with 
^ = 1.65 X 10^^~^ and orange with ^ = 4.85 x 10^ to show the evolution of the 
inversion within the cavity as a function of the pump strength. 



D. Simulation constants 



In this appendix we list the parameters used in each of the FDTD simulations described above. 
These constants are matrices, with Yij denoting the decay rate from |j) to |/). These values are 
given in their dimensionless form, i.e. ymeas = 7rQa\L/c. Unlisted entries are zero. We also note 
that throughout this paper |0) denotes the ground state, so these matrices are indexed. 
For the four-level simulations in Fig. [51 



741v, fig|2] = 



0.8 



5x 10" 
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(47) 



The dipole matrix element is g = 2.3 • 10 



-i2j^3/2^ and the number of gain atoms is n = 5 - 
j^Q23j^-3 jYiQ pump ^ was varied between 3 • 10~^ and 3 • 10~^. 

Thus, for an optical wavelength of A = 628nm, the requirement in Fig. [2l that nokL = 60 
means that L = 4/im. Using this length, the decay rates can be converted to their unit- full values 

as Y± 
^ = 3- 



3 • lO^^s"^ 723 = 701 = 6 • lO^-^s-^ 7i2 = 3.75 • lO^^s"^ and the pump at threshold is 
10^s~^ Similarly, the dipole matrix element also acquires units of inverse time, and 
can be expressed as g'^/h = 3.98 • 10~^m^/s, which corresponds to a coupling constant in the 
classical oscillating dipole picture of a = 10~^C^/kg. These constants can be seen to be similar 
to those used in other studies of optical microcavities ll4lfT2l. 



For the six-level simulations in Fig. [21 
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Furthermore, 715 = 10-^, and the lasing transition is between levels |3) and |1) (where the 
ground state is again |0) and the states are numbered in order of increasing energy). 
For the four-level simulations in Fig.O 
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For the six-level simulations in Fig.[3l 
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The four-level simulations of the random cavity in Fig. |4] used the same parameters as the 
four-level simulations in Fig.[2l 



